
# tif files are downloaded from https://www.ngdc.noaa.gov/eog/data/web_data/v4composites/
# extracted from tar and put in directory data/tiff/
# run to crate input data for read_lights_drought.do

#first set working directory on level down from data
#setwd()

# Start the clock
ptm <- proc.time()

library(raster)


for (i in 2000:2011) {


infile <- paste("./data/tiff/F", i, ".v4b_web.stable_lights.avg_vis.tif", sep="")
outfile <- paste("./data/grid_ascii/lights_05deg_", i, ".asc", sep="")
outfile_ihs <- paste("./data/grid_ascii/lights_05deg_ihs_", i, ".asc", sep="")

r <- raster(infile)

#levels
#aggregate to .05 degrees
ra <- aggregate(r, fact=6, na.rm=TRUE, fun=mean)

#crop at lat=-60 to match droughts
rc <- crop(ra, extent(ra, 1, 2700, 1, 7200))

writeRaster(rc, filename=outfile, format = "ascii", overwrite=TRUE, NAflag=-9999)


#inverse hyperbolic sine
#aggregate to .05 degrees
ra_ihs <- aggregate(asinh(r), fact=6, na.rm=TRUE, fun=mean)

#crop at lat=-60 to match droughts
rc_ihs <- crop(ra_ihs, extent(ra, 1, 2700, 1, 7200))

writeRaster(rc_ihs, filename=outfile_ihs, format = "ascii", overwrite=TRUE, NAflag=-9999)
}

# Stop the clock
proc.time() - ptm

